######################################
######################################
######################################
###Coefficient Impact Plots Percent###
######################################
######################################
######################################
library(mvtnorm)
#####################
#####################
###State Conflicts###
#####################
#####################

##############
###Baseline###
##############
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","p_anom","spei","outbreak","t_anom","gid","month")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
p1<-NULL
p2<-NULL
all.expectedvalues1<-NULL
all.expectedvalues2<-NULL
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm2.state))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm2.state))
  xoutcome.1<-rbind(log.nl,logged.pop,p_anom,spei,0,t_anom)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,p_anom,spei,1,t_anom)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.st2 <- history
hist.per.st2 <- hist.per
history.st2
hist.per.st2


##########
###Full###
##########
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","lagacled_battle_state","t","p_anom","spei","outbreak","t_anom","gid","month","logLifeExpectancy","GovEff","logGDPpc")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
lagacled_battle_state=mean(full.na$lagacled_battle_state)
t=mean(full.na$t)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
logLifeExpectancy=mean(full.na$logLifeExpectancy)
GovEff=mean(full.na$GovEff)
logGDPpc=mean(full.na$logGDPpc)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm3.state))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm3.state))
  xoutcome.1<-rbind(log.nl,logged.pop,lagacled_battle_state,t,p_anom,spei,0,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,lagacled_battle_state,t,p_anom,spei,1,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.st3 <- history
hist.per.st3 <- hist.per
history.st3
hist.per.st3


#####################
#####################
###Rebel Conflicts###
#####################
#####################

##############
###Baseline###
##############
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","p_anom","spei","outbreak","t_anom","gid","month")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
p1<-NULL
p2<-NULL
all.expectedvalues1<-NULL
all.expectedvalues2<-NULL
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm2.rebel))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm2.rebel))
  xoutcome.1<-rbind(log.nl,logged.pop,p_anom,spei,0,t_anom)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,p_anom,spei,1,t_anom)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.reb2 <- history
hist.per.reb2 <- hist.per
history.reb2
hist.per.reb2


##########
###Full###
##########
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","lagacled_battle_rebel","t","p_anom","spei","outbreak","t_anom","gid","month","logLifeExpectancy","GovEff","logGDPpc")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
lagacled_battle_rebel=mean(full.na$lagacled_battle_rebel)
t=mean(full.na$t)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
logLifeExpectancy=mean(full.na$logLifeExpectancy)
GovEff=mean(full.na$GovEff)
logGDPpc=mean(full.na$logGDPpc)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm3.rebel))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm3.rebel))
  xoutcome.1<-rbind(log.nl,logged.pop,lagacled_battle_rebel,t,p_anom,spei,0,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,lagacled_battle_rebel,t,p_anom,spei,1,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.reb3 <- history
hist.per.reb3 <- hist.per
history.reb3
hist.per.reb3


#################################
#################################
###Political Militia Conflicts###
#################################
#################################

##############
###Baseline###
##############
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","p_anom","spei","outbreak","t_anom","gid","month")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
p1<-NULL
p2<-NULL
all.expectedvalues1<-NULL
all.expectedvalues2<-NULL
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm2.mil.pol))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm2.mil.pol))
  xoutcome.1<-rbind(log.nl,logged.pop,p_anom,spei,0,t_anom)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,p_anom,spei,1,t_anom)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.polmil2 <- history
hist.per.polmil2 <- hist.per
history.polmil2
hist.per.polmil2


##########
###Full###
##########
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","lagacled_battle_polmil","t","p_anom","spei","outbreak","t_anom","gid","month","logLifeExpectancy","GovEff","logGDPpc")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
lagacled_battle_polmil=mean(full.na$lagacled_battle_polmil)
t=mean(full.na$t)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
logLifeExpectancy=mean(full.na$logLifeExpectancy)
GovEff=mean(full.na$GovEff)
logGDPpc=mean(full.na$logGDPpc)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm3.mil.pol))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm3.mil.pol))
  xoutcome.1<-rbind(log.nl,logged.pop,lagacled_battle_polmil,t,p_anom,spei,0,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,lagacled_battle_polmil,t,p_anom,spei,1,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.polmil3 <- history
hist.per.polmil3 <- hist.per
history.polmil3
hist.per.polmil3


#################################
#################################
###Identity Militia Conflicts###
#################################
#################################

##############
###Baseline###
##############
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","p_anom","spei","outbreak","t_anom","gid","month")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
p1<-NULL
p2<-NULL
all.expectedvalues1<-NULL
all.expectedvalues2<-NULL
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm2.mil.id))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm2.mil.id))
  xoutcome.1<-rbind(log.nl,logged.pop,p_anom,spei,0,t_anom)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,p_anom,spei,1,t_anom)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.idmil2 <- history
hist.per.idmil2 <- hist.per
history.idmil2
hist.per.idmil2


##########
###Full###
##########
#Subset NA data
full.na <- na.omit(full.1996.2019[,c("log.nl","logged.pop","lagacled_battle_idmil","t","p_anom","spei","outbreak","t_anom","gid","month","logLifeExpectancy","GovEff","logGDPpc")])
###Create indicator terms
log.nl=mean(full.na$log.nl)
logged.pop=mean(full.na$logged.pop)
lagacled_battle_idmil=mean(full.na$lagacled_battle_idmil)
t=mean(full.na$t)
p_anom=mean(full.na$p_anom)
spei=mean(full.na$spei)
#out
t_anom=mean(full.na$t_anom)
logLifeExpectancy=mean(full.na$logLifeExpectancy)
GovEff=mean(full.na$GovEff)
logGDPpc=mean(full.na$logGDPpc)
###Set seed
set.seed(1234)
#set sims  (number of simulations)
sims <- 10000
#create some storage matrices
First.Differences<-NULL
First.Diff.Per <- NULL
all.first.diff.per <- NULL
all.firstdifferences<-NULL
history<-matrix(0,nrow=1,ncol=3)
hist.per <- matrix(0,nrow=1,ncol=3)
#Set coefficients
c <- c(coefficients(lm3.mil.id))
#start loop
for(j in 1:sims){ 
  #predicted#
  bsim<-rmvnorm(1, mean=c, sigma=vcov(lm3.mil.id))
  xoutcome.1<-rbind(log.nl,logged.pop,lagacled_battle_idmil,t,p_anom,spei,0,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ####Creating the matrices
  equationcount1<-bsim %*% xoutcome.1
  ###First difference change prediction
  xoutcome.2<-rbind(log.nl,logged.pop,lagacled_battle_idmil,t,p_anom,spei,1,t_anom,
                    logLifeExpectancy,GovEff,logGDPpc)
  ###Creating the matrices
  equationcount2<-bsim %*% xoutcome.2
  Firstdifference<-equationcount2-equationcount1       
  First.Diff.Per <- ((equationcount2-equationcount1)/abs(equationcount1))*100
  all.firstdifferences<-c(all.firstdifferences, Firstdifference)
  all.first.diff.per<- c(all.first.diff.per, First.Diff.Per)
  ###end loop
}
###Create summary
sorted.firstdifferences<-sort(all.firstdifferences)
history[1,1]<-mean(sorted.firstdifferences)
history[1,2]<-sorted.firstdifferences[round(0.95*length(sorted.firstdifferences))]
history[1,3]<-sorted.firstdifferences[round(0.05*length(sorted.firstdifferences))]
###Percent summary
sorted.first.diff.per<-sort(all.first.diff.per)
hist.per[1,1]<- mean(sorted.first.diff.per)
hist.per[1,2] <- sorted.first.diff.per[round(0.95*length(sorted.first.diff.per))]
hist.per[1,3] <- sorted.first.diff.per[round(0.05*length(sorted.first.diff.per))]
###end final loop
#}
history.idmil3 <- history
hist.per.idmil3 <- hist.per
history.idmil3
hist.per.idmil3



#####################
#####################
###Generate Plots###
#####################
#####################


##########################
###Baseline Model Plots###
##########################

#Mean
plot.mean <- c(hist.per.st2[,1],hist.per.reb2[,1],hist.per.polmil2[,1],hist.per.idmil2[,1])
#Upper bound
plot.up <- c(hist.per.st2[,2],hist.per.reb2[,2],hist.per.polmil2[,2],hist.per.idmil2[,2])
#Lower bound
plot.low <- c(hist.per.st2[,3],hist.per.reb2[,3],hist.per.polmil2[,3],hist.per.idmil2[,3])

##########
###Plot###

###Set up plotting space
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-210,210)
xlim<-c(0,5)
y<-plot.mean
x<-seq(1:length(y))
yupper<-plot.up
ylower<-plot.low


###plot actual effects
setwd("~/OneDrive - Indiana University/FromGoogle/GeolocatedDisease/Revision/")
#open file
jpeg("fig1.jpeg", width = 13.59, height = 10.02, units = 'in', res = 500)
par(mfrow=c(2,1))
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col=ifelse((x==1|x==4), "black", "darkgrey"),pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col=ifelse((x==1|x==4), "black", "darkgrey"),lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=1)

#Add sectioning lines
abline(v=2.5,lty=3)

#Label X-axis values
axis(1, at=c(1,2,3,4), lab=c("State","Rebel","Pol. Mil.", "Id. Mil."), tick = F, cex.axis=0.9, font=2)

axis(3, at=c(1.5,3.5), lab=c("Armed conflict","Social conflict"), tick = TRUE, cex.axis=1.1, font=2)

#axis(1, at=c(1,2,3,4,5,6,7,8,9,10), lab=c("(SB)","(SF)","(RB)","(RF)","(MB)","(MF)","(PMB)","(PMF)", "(IMB)", "(IMF)"), tick = TRUE, cex.axis=0.9)

#Add a title
title(ylab="% Change in Expected Conflict Rates", main="") #, cex.main=1)

#Adjust Y-Axis Labels
axis(2, at=c((-200), (0), 200, 400), lab=c("-200%",  "0", "200%", "400%"), las=2, cex.axis=0.8)

#Add box around plot
box()


######################
###Full Model Plots###
######################

#Mean
plot.mean <- c(hist.per.st3[,1],hist.per.reb3[,1],hist.per.polmil3[,1],hist.per.idmil3[,1])
#Upper bound
plot.up <- c(hist.per.st3[,2],hist.per.reb3[,2],hist.per.polmil3[,2],hist.per.idmil3[,2])
#Lower bound
plot.low <- c(hist.per.st3[,3],hist.per.reb3[,3],hist.per.polmil3[,3],hist.per.idmil3[,3])

##########
###Plot###

###Set up plotting space
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-30,140)
xlim<-c(0,5)
y<-plot.mean
x<-seq(1:length(y))
yupper<-plot.up
ylower<-plot.low


###plot actual effects
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col=ifelse((x==1|x==4), "black", "darkgrey"),pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col=ifelse((x==1|x==4), "black", "darkgrey"),lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=1)

#Add sectioning lines
#lines(x, ave.y, col="black", lty="2")
abline(v=2.5,lty=3)
#abline(v=4.5,lty=3)
#abline(v=6.5,lty=3)

#Label X-axis values
axis(1, at=c(1,2,3,4), lab=c("State","Rebel","Pol. Mil.", "Id. Mil."), tick = F, cex.axis=0.9, font=2)

axis(3, at=c(1.5,3.5), lab=c("Armed conflict","Social conflict"), tick = TRUE, cex.axis=1.1, font=2)

#Add a title
title(ylab="% Change in Expected Conflict Rates", main="") #, cex.main=1)

#Adjust Y-Axis Labels
axis(2, at=c((-30), (0), 60, 120), lab=c("-30%",  "0", "60%", "120%"), las=2, cex.axis=0.8)

#Add box around plot
box()
#close file
dev.off()

